#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/times.h>
#include <omp.h>
#include "mytime.h"

#define N_DEF 200
#define TOL_DEF 0.0000001
#define MAXIT_DEF 500

/****************/
int *despl;
int *count;
/****************/

double
jacobi_omp (double **a, int n, double *b, double *x, double tol, int maxit)
{
  int i, j, it = 0;
  double sum, resi = 10.0, *new_x;

  new_x = (double *) malloc (n * sizeof (double));

  for (i = 0; i < n; i++)
    x[i] = b[i];
  
  while (it < maxit)
  {
    resi = 0.0;
/* #pragma omp target device(mpi) */
#pragma omp parallel for private (sum, j) reduction (+ : resi)
/* #pragma llc result (&new_x[i], 1) */
    for (i = 0; i < n; i++)
	{
	  sum = 0.0;
	  for (j = 0; j < n; j++)
	    sum += a[i][j] * x[j];
	  resi += fabs (sum - b[i]);
	  sum += -a[i][i] * x[i];
	  new_x[i] = (b[i] - sum) / a[i][i];
	}
    for (i = 0; i < n; i++)
	  x[i] = new_x[i];
    it++;
  }

  free (new_x);

  return resi;
}

int
main (int argc, char *argv[])
{
  double tol;
  double **a, *b, *x;
  int i, j, n; 
  // int nfilas;
  int diag, maxit;
  int nprocs, miid;

  CLOCK_TYPE chrono;
  double omp_time;
  double omp_resi;

  n = N_DEF;
  maxit = MAXIT_DEF;
  tol = TOL_DEF;

  switch (argc)
  {
    case 4:
      tol = atof (argv[3]);
    case 3:
      maxit = atoi (argv[2]);
    case 2:
      n = atoi (argv[1]);
  }


  /* Por simplicidad supondremos que n es divisible por el numero de procesos: nprocs */
  // nfilas = n / nprocs;

  ReservaMatriz (&a, n, n);
  ReservaVector (&b, n);
  ReservaVector (&x, n);

#pragma omp master
  {
    GeneraVector (b, n);
    GeneraMatriz (a, n, n);
    /* Para garantizar la convergencia del metodo de Jacobi */
    /* Convierte en diagonal dominante la matriz distribuida */
    diag = 0;
    for (i = 0; i < n; i++)
    {
	  for (j = 0; j < n; j++)
	    a[i][diag] += a[i][j];
	  diag++;
    }
  }
  /* envio de la matriz */
  //for (i = 0; i < n; i++)
  //  MPI_Bcast (a[i], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  //MPI_Barrier (MPI_COMM_WORLD);
  #pragma omp barrier

  /****************************************************************************
   jacobi omp */
  CLOCK_Start (chrono);
  omp_resi = jacobi_omp (a, n, b, x, tol, maxit);
  CLOCK_End (chrono, omp_time);
  //MPI_Barrier (MPI_COMM_WORLD);
  #pragma omp barrier

/*  CompruebaSolucion (a, x, b, n);*/

  /****************************************************************************
  ****************************************************************************/
 
  #pragma omp parallel private(miid, nprocs)
  {
    nprocs = omp_get_num_threads();
    miid = omp_get_thread_num();
    printf ("%d:%g:%g\n", miid, omp_time, omp_resi);
  }

  /**************************************************************************
   Libeacion de memoria */
  LiberaMatriz (a);
  LiberaVector (b);
  LiberaVector (x);

  return 0;

}
